##########################################################
### Disorganized Political Violence:                   ###
### A Demonstration Case of Temperature and Insurgency ###
### Replication Code: Primary Results --               ###
### Temperature Deviations Results, w/                 ###
### GAM and Binned Results (Appendix Results)          ###
### Andrew Shaver, Alex Bollfrass                      ###
### Date: 07/04/2022                                   ###
##########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(MASS)
library(lubridate)
library(DataCombine)
library(broom)
library(data.table)
library(radiant.model)
library(lfe)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(lfe)
library(arm)
library(mgcv)
library(latex2exp)

packageVersion("MASS")
# ‘7.3.55’
packageVersion("lubridate")
# ‘1.8.0’
packageVersion("DataCombine")
# ‘0.2.21’
packageVersion("broom")
# ‘0.7.12’
packageVersion("data.table")
# ‘1.14.2’
packageVersion("radiant.model")
# ‘1.4.4’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("arm")
# ‘1.12.2’
packageVersion("mgcv")
# ‘1.8.39’
packageVersion("latex2exp")
# ‘0.9.4’

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load data:
# For primary temperature-violence analysis:
# Daily analysis:
iq_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iraq_baghdad_basrah_processed_data.rds")
af_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/afghanstan_viol_districts_processed_data.rds")
# For attitude analysis:
bg_sur <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/baghdad_survey.rds")
# For median analysis:
iq_median_temp <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/Affective_Motivations/MajorRevision/R3_Submission/ProcessedData/iq_daily_median.rds")

# Match median temperature data to the Iraq panel:
colnames(iq_median_temp)[2] <- "governorate"
iq_sig1 <- plyr::join(iq_sig1, iq_median_temp, by = c("date", "governorate"))

# Next, construct temperature deviations measure:
# Write function:
deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(max ~ lp(nday, nn = .05), data=temp)
    temp$temp_hat <- predict(t, newdata = temp$nday)
    temp$temp_dev <- temp$max - temp$temp_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_sig1
iq_sig1 <- deviation("governorate")

panel <- af_sig1
af_sig1 <- deviation("station")

iq_sig1 <- as.data.frame(iq_sig1)
af_sig1 <- as.data.frame(af_sig1)

deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(median ~ lp(nday, nn = .05), data=temp)
    temp$median_hat <- predict(t, newdata = temp$nday)
    temp$median_dev <- temp$median - temp$median_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_sig1
iq_sig1 <- deviation("governorate")
iq_sig1 <- as.data.frame(iq_sig1)

# For temperature deviations for the Baghdad survey data, given that the unit is not the
# match deviations from panel:
# Columns to keep:
to_keep <- c("date", "temp_dev", "temp_1", "temp_2", "temp_3", "temp_4", 
             "temp_5", "temp_6", "temp_7", "daylight")
bg_sur <- plyr::join(bg_sur, iq_sig1[iq_sig1$governorate %in% "Baghdad", to_keep], by = "date")

# Create survey dataset removing top income earners:
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]

# Next, generate descriptive statistics (TABLES 1, 2 AND 3 IN PAPER):
stargazer(iq_sig1[, c("governorate", "date", "max", "prcp", "wdsp", "dewp", "visib", "mxspd", "daylight", "hoursofpower", 
                      "unconstrained", "idf", "ied", "constrained")], iqr = T)
stargazer(af_sig1[, c("district", "date", "max", "prcp", "wdsp", "dewp", "visib", "mxspd", "daylight", 
                      "df", "idf", "ied")], iqr = T)
stargazer(bg_sur[, c("block", "date", "max", "prcp", "wdsp", "dewp", "visib", "mxspd", "daylight", "age",
                     "work_wkhr", "income.w", "hhld_size", "m_sig1_7", "p3_electr_binary")], iqr = T)


# Next, generate temperature deviations for each context:
# Set temperature range for all:
temperature_range <- 75:105

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# OLS model results:

# Check the normality of the errors:
iq_sig1$threshold <- ifelse(iq_sig1$max > 75, 1, 0)
test1 <- lm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(governorate), data = iq_sig1)
test2 <- lm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(governorate), data = iq_sig1)

test1.stdres = rstandard(test1)
qqnorm(test1.stdres)
qqline(test1.stdres)

test2.stdres = rstandard(test2)
qqnorm(test2.stdres)
qqline(test2.stdres)

temperature_deviation_iq_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols <- temperature_deviation_iq_ols(iq_sig1, 0.95)
iq_90_ols <- temperature_deviation_iq_ols(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols$estimate) / mean(log(iq_sig1$unconstrained+1))

# For plotting Baghdad-Basrah results:
temp_ranges <- unique(iq_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

# Plot:
p_iq_ols <- ggplot() + 
  
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 75:105, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Negative binomial model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_iq <- data.frame(temp_dev = c(0, sd(iq_sig1$temp_dev)),
                          threshold = 0,
                          
                          prcp = mean(iq_sig1$prcp, na.rm = T), 
                          wdsp = mean(iq_sig1$wdsp, na.rm = T), 
                          dewp = mean(iq_sig1$dewp, na.rm = T), 
                          visib = mean(iq_sig1$visib, na.rm = T), 
                          mxspd = mean(iq_sig1$mxspd, na.rm = T), 
                          daylight = mean(iq_sig1$daylight, na.rm = T), 
                          hoursofpower = mean(iq_sig1$hoursofpower, na.rm = T), 
                          
                          idf = mean(iq_sig1$idf), 
                          ied = mean(iq_sig1$ied),  
                          constrained = mean(iq_sig1$constrained),  
                          unconstrained = mean(iq_sig1$unconstrained),  
                          
                          temp_1 = mean(iq_sig1$temp_1, na.rm = T),
                          temp_2 = mean(iq_sig1$temp_2, na.rm = T),
                          temp_3 = mean(iq_sig1$temp_3, na.rm = T),
                          temp_4 = mean(iq_sig1$temp_4, na.rm = T),
                          temp_5 = mean(iq_sig1$temp_5, na.rm = T),
                          temp_6 = mean(iq_sig1$temp_6, na.rm = T),
                          temp_7 = mean(iq_sig1$temp_7, na.rm = T),
                          
                          unconstrained_1 = mean(iq_sig1$unconstrained_1, na.rm = T),
                          unconstrained_2 = mean(iq_sig1$unconstrained_2, na.rm = T),
                          unconstrained_3 = mean(iq_sig1$unconstrained_3, na.rm = T),
                          unconstrained_4 = mean(iq_sig1$unconstrained_4, na.rm = T),
                          unconstrained_5 = mean(iq_sig1$unconstrained_5, na.rm = T),
                          unconstrained_6 = mean(iq_sig1$unconstrained_6, na.rm = T),
                          unconstrained_7 = mean(iq_sig1$unconstrained_7, na.rm = T),
                          
                          constrained_1 = mean(iq_sig1$constrained_1, na.rm = T),
                          constrained_2 = mean(iq_sig1$constrained_2, na.rm = T),
                          constrained_3 = mean(iq_sig1$constrained_3, na.rm = T),
                          constrained_4 = mean(iq_sig1$constrained_4, na.rm = T),
                          constrained_5 = mean(iq_sig1$constrained_5, na.rm = T),
                          constrained_6 = mean(iq_sig1$constrained_6, na.rm = T),
                          constrained_7 = mean(iq_sig1$constrained_7, na.rm = T),
                          
                          idf_1 = mean(iq_sig1$idf_1, na.rm = T),
                          idf_2 = mean(iq_sig1$idf_2, na.rm = T),
                          idf_3 = mean(iq_sig1$idf_3, na.rm = T),
                          idf_4 = mean(iq_sig1$idf_4, na.rm = T),
                          idf_5 = mean(iq_sig1$idf_5, na.rm = T),
                          idf_6 = mean(iq_sig1$idf_6, na.rm = T),
                          idf_7 = mean(iq_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_sig1$ied_1, na.rm = T),
                          ied_2 = mean(iq_sig1$ied_2, na.rm = T),
                          ied_3 = mean(iq_sig1$ied_3, na.rm = T),
                          ied_4 = mean(iq_sig1$ied_4, na.rm = T),
                          ied_5 = mean(iq_sig1$ied_5, na.rm = T),
                          ied_6 = mean(iq_sig1$ied_6, na.rm = T),
                          ied_7 = mean(iq_sig1$ied_7, na.rm = T),
                          
                          week = median(iq_sig1$week),
                          month = median(iq_sig1$month),
                          governorate = "Baghdad")

temperature_deviation_iq_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$unconstrained)
    
    reg_temp_month.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb <- temperature_deviation_iq_nb(iq_sig1, 0.95)
iq_90_nb <- temperature_deviation_iq_nb(iq_sig1, 0.90)

# Plot:
temp <- scales::percent(iq_95_nb[iq_95_nb$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_iq_nb <- ggplot() + 
  
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb[iq_95_nb$time %in% "week", ]$gamma, labels=temp)) + 

  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 75:105, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.025, vjust=-6.5, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y= 0.045, label="Unconstrained",
           color="black", size = 6, fontface = 2, hjust = 0)

# Constrained Violence:
# OLS model results:
temperature_deviation_iq_ols_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.95)
iq_90_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_con$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols_con$estimate) / mean(log(iq_sig1$constrained+1))

# Plot:
p_iq_ols_con <- ggplot() + 
  
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "green4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 75:105, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$constrained)
    
    reg_temp_month.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.95)
iq_90_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.90)

# Plot:
temp <- scales::percent(iq_95_nb_con[iq_95_nb_con$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_iq_nb_con <- ggplot() + 
  
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "green4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb_con[iq_95_nb_con$time %in% "week", ]$gamma, labels=temp)) + 

  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 75:105, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.025, vjust=-7.25, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y= 0.045, label="Constrained",
           color="black", size = 6, fontface = 2, hjust = 0)

# IED Attacks:
# OLS model results:
temperature_deviation_iq_ols_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.95)
iq_90_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_ied$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols_ied$estimate) / mean(log(iq_sig1$ied+1))

# Plot:
p_iq_ols_ied <- ggplot() + 
  
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 75:105, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$ied)
    
    reg_temp_month.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.95)
iq_90_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.90)

# Plot:
temp <- scales::percent(iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_iq_nb_ied <- ggplot() + 
  
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ]$gamma, labels=temp)) + 

  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 75:105, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.0275, vjust=-6.75, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y= 0.045, label="IED",
           color="black", size = 6, fontface = 2, hjust = 0) + 
  annotate("point", x = 89, y = 0.0425, colour = "blue4", size = 3) + 
  annotate("point", x = 90, y = 0.0425, colour = "green4", size = 3) + 
  annotate("point", x = 91, y = 0.0425, colour = "firebrick", size = 3) + 
  
  annotate(geom="text", x=97, y= 0.0425, label="Primary Specifications",
           color="black", size = 3) + 
  annotate("point", x = 89, y = 0.035, colour = "gray60", size = 3) + 
  annotate(geom="text", x=96.5, y= 0.035, label="Robustness Checks",
           color="black", size = 3) +
  geom_rect(aes(xmin = 87.5, xmax = 102.5, ymin = 0.0275, ymax = 0.0475),
            fill = "transparent", color = "black", size = 0.5)

iraq_combined <- ggarrange(p_iq_nb + xlab("") + ylab("Quasi-Poisson (Count) Models"), 
                           p_iq_nb_con + xlab("") + ylab(""), 
                           p_iq_nb_ied + xlab("") + ylab(""),  
                           p_iq_ols + xlab("") + ylab("OLS Regression Models"),
                           p_iq_ols_con + xlab("") + ylab(""), 
                           p_iq_ols_ied + xlab("") + ylab(""), 
                           ncol = 3, nrow=2, labels  = "")

# PRODUCES FIGURE 21 IN PAPER:
annotate_figure(iraq_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Violence Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

iraq_combined <- ggarrange(p_iq_nb + xlab("") + ylab("Quasi-Poisson (Count) Models"), 
                           p_iq_nb_con + xlab("") + ylab(""), 
                           p_iq_nb_ied + xlab("") + ylab(""),  
                           ncol = 3, nrow=1, labels  = c("", "", ""))

# PRODUCES FIGURE 3 IN PAPER:
annotate_figure(iraq_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Violence Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### AFGHANISTAN.
###

# Afghanistan districts:
# OLS:

# Check the normality of the errors:
af_sig1$threshold <- ifelse(af_sig1$max > 75, 1, 0)
test1 <- lm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(station), data = af_sig1)
test2 <- lm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(station), data = af_sig1)

test1.stdres = rstandard(test1)
qqnorm(test1.stdres)
qqline(test1.stdres)

test2.stdres = rstandard(test2)
qqnorm(test2.stdres)
qqline(test2.stdres)

# DF:
# Log:
temperature_deviation_af_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + station | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + station | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                                 week + station | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                                  month + station | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_ols <- temperature_deviation_af_ols(af_sig1, 0.95)
af_90_ols <- temperature_deviation_af_ols(af_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95_ols$change <- (sd(af_sig1$temp_dev, na.rm = T) * af_95_ols$estimate) / mean(log(af_sig1$df+1))

# For plotting Afghan district results:
temp_ranges_af <- unique(af_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges_af)*2-2), by=2)){
  temp_ranges_af <- append(temp_ranges_af, " ", after=i)
}

p_af_ols <- ggplot() + 
  
  geom_line(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 75:105, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Negative binomial model:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_af <- data.frame(temp_dev = c(0, sd(af_sig1$temp_dev)),
                          threshold = 0,
                          
                          prcp = mean(af_sig1$prcp, na.rm = T), 
                          wdsp = mean(af_sig1$wdsp, na.rm = T), 
                          dewp = mean(af_sig1$dewp, na.rm = T), 
                          visib = mean(af_sig1$visib, na.rm = T), 
                          mxspd = mean(af_sig1$mxspd, na.rm = T), 
                          daylight = mean(af_sig1$daylight, na.rm = T), 
                          
                          idf = mean(af_sig1$idf), 
                          ied = mean(af_sig1$ied),  
                          df = mean(af_sig1$df), 
                          
                          temp_1 = mean(af_sig1$temp_1, na.rm = T),
                          temp_2 = mean(af_sig1$temp_2, na.rm = T),
                          temp_3 = mean(af_sig1$temp_3, na.rm = T),
                          temp_4 = mean(af_sig1$temp_4, na.rm = T),
                          temp_5 = mean(af_sig1$temp_5, na.rm = T),
                          temp_6 = mean(af_sig1$temp_6, na.rm = T),
                          temp_7 = mean(af_sig1$temp_7, na.rm = T),
                          
                          df_1 = mean(af_sig1$df_1, na.rm = T),
                          df_2 = mean(af_sig1$df_2, na.rm = T),
                          df_3 = mean(af_sig1$df_3, na.rm = T),
                          df_4 = mean(af_sig1$df_4, na.rm = T),
                          df_5 = mean(af_sig1$df_5, na.rm = T),
                          df_6 = mean(af_sig1$df_6, na.rm = T),
                          df_7 = mean(af_sig1$df_7, na.rm = T),
                          
                          idf_1 = mean(af_sig1$idf_1, na.rm = T),
                          idf_2 = mean(af_sig1$idf_2, na.rm = T),
                          idf_3 = mean(af_sig1$idf_3, na.rm = T),
                          idf_4 = mean(af_sig1$idf_4, na.rm = T),
                          idf_5 = mean(af_sig1$idf_5, na.rm = T),
                          idf_6 = mean(af_sig1$idf_6, na.rm = T),
                          idf_7 = mean(af_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(af_sig1$ied_1, na.rm = T),
                          ied_2 = mean(af_sig1$ied_2, na.rm = T),
                          ied_3 = mean(af_sig1$ied_3, na.rm = T),
                          ied_4 = mean(af_sig1$ied_4, na.rm = T),
                          ied_5 = mean(af_sig1$ied_5, na.rm = T),
                          ied_6 = mean(af_sig1$ied_6, na.rm = T),
                          ied_7 = mean(af_sig1$ied_7, na.rm = T),
                          
                          week = median(af_sig1$week),
                          month = median(af_sig1$month),
                          station = "409900")

temperature_deviation_af_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_af[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_af[1,], type = "response")) /  mean(af_sig1$df)
    
    reg_temp_month.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                                   as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                                    as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_nb <- temperature_deviation_af_nb(af_sig1, 0.95)
af_90_nb <- temperature_deviation_af_nb(af_sig1, 0.90)

temp <- scales::percent(af_95_nb[af_95_nb$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_af_qp <- ggplot() + 
  
  geom_line(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = af_95_nb[af_95_nb$time %in% "week", ]$gamma, labels=temp)) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 75:105, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.025, vjust=-6.5, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y=0.016, label="Direct Fire",
           color="black", size = 6, fontface = 2, hjust = 0)

# IED:
temperature_deviation_af_ols_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + station | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + station | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df |
                                 week + station | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df |
                                  month + station | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_ols_ied <- temperature_deviation_af_ols_ied(af_sig1, 0.95)
af_90_ols_ied <- temperature_deviation_af_ols_ied(af_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95_ols_ied$change <- (sd(af_sig1$temp_dev, na.rm = T) * af_95_ols_ied$estimate) / mean(log(af_sig1$ied+1))

p_af_ols_ied <- ggplot() + 
  
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 75:105, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# QP:
temperature_deviation_af_nb_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_af[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_af[1,], type = "response")) /  mean(af_sig1$df)
    
    reg_temp_month.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                                   as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                                    as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_nb_ied <- temperature_deviation_af_nb_ied(af_sig1, 0.95)
af_90_nb_ied <- temperature_deviation_af_nb_ied(af_sig1, 0.90)

temp <- scales::percent(af_95_nb_ied[af_95_nb_ied$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_af_qp_ied <- ggplot() + 
  
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = af_95_nb_ied[iq_95_nb_ied$time %in% "week", ]$gamma, labels=temp)) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 75:105, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.0275, vjust=-6.75, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y=0.016, label="IED",
           color="black", size = 6, fontface = 2, hjust = 0) + 
  annotate("point", x = 91, y = 0.015, colour = "blue4", size = 3) + 
  annotate("point", x = 92, y = 0.015, colour = "firebrick", size = 3) + 
  
  annotate(geom="text", x=97, y= 0.015, label="Primary Specifications",
           color="black", size = 3) + 
  annotate("point", x = 91, y = 0.013, colour = "gray60", size = 3) + 
  annotate(geom="text", x=96.6, y= 0.013, label="Robustness Checks",
           color="black", size = 3) +
  geom_rect(aes(xmin = 90, xmax = 101.5, ymin = 0.0115, ymax = 0.0165),
            fill = "transparent", color = "black", size = 0.5)

afghanistan_combined <- ggarrange(p_af_qp + xlab("") + ylab("Quasi-Poisson (Count) Models"),
                                  p_af_qp_ied + xlab("") + ylab(""),
                                  p_af_ols + xlab("") + ylab("OLS Regression Models"),
                                  p_af_ols_ied + xlab("") + ylab(""),
                                  ncol = 2, nrow=2, labels = "")

# PRODUCES FIGURE 22 IN PAPER:
annotate_figure(afghanistan_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Violence Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

afghanistan_combined <- ggarrange(p_af_qp + xlab("") + ylab("Quasi-Poisson (Count) Models"),
                                  p_af_qp_ied + xlab("") + ylab(""),
                                  ncol = 2, nrow=1, labels = c("", ""))

# PRODUCES FIGURE 4 IN PAPER:
annotate_figure(afghanistan_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Violence Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### BAGHDAD SURVEYS.
###

# OLS model results:
temperature_deviation_bg_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                            p3_electr_binary +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                            education + mahala + date_month | 0 | mahala, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                             p3_electr_binary + 
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                             education + block + date_month | 0 | block, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight |
                                 p3_electr_binary +
                                 education + mahala + date_month | 0 | mahala, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                                  p3_electr_binary |
                                  education + block + date_month | 0 | block, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

bg_95_ols <- temperature_deviation_bg_ols(bg_sur.s3.l, 0.95)
bg_90_ols <- temperature_deviation_bg_ols(bg_sur.s3.l, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
bg_95_ols$change <- (sd(bg_sur.s3.l$temp_dev, na.rm = T) * bg_95_ols$estimate) / mean(bg_sur.s3.l$vs_mnf)

# For plotting Baghdad survey results:
temp_ranges_bg <- unique(bg_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges_bg)*2-2), by=2)){
  temp_ranges_bg <- append(temp_ranges_bg, " ", after=i)
}

p_bg_ols <- ggplot() +
  
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size =  0.5, alpha = .35, color = "gray60") +
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 1, color = "gray60") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  ylab(expression(beta~Delta*"Temperature")) +
  xlab(expression(gamma*" (°F)")) +
  ylim(-0.004, 0.0125) +
  annotate(geom = "text", x = 75:105, y = -0.00325,
           label = temp_ranges_bg,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") +
  theme_bw() +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Logit model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_bg <- data.frame(temp_dev = c(0, sd(bg_sur.s3.l$temp_dev)),
                          threshold = 0,
                          
                          prcp = mean(bg_sur.s3.l$prcp, na.rm = T), 
                          wdsp = mean(bg_sur.s3.l$wdsp, na.rm = T), 
                          dewp = mean(bg_sur.s3.l$dewp, na.rm = T), 
                          visib = mean(bg_sur.s3.l$visib, na.rm = T), 
                          mxspd = mean(bg_sur.s3.l$mxspd, na.rm = T), 
                          daylight = mean(bg_sur.s3.l$daylight, na.rm = T), 
                          
                          age = mean(bg_sur.s3.l$age, na.rm = T),  
                          work_wkhr = mean(bg_sur.s3.l$work_wkhr),  
                          income.w = mean(bg_sur.s3.l$income.w),  
                          hhld_size = mean(bg_sur.s3.l$hhld_size), 
                          m_sig1_7 = mean(bg_sur.s3.l$m_sig1_7),  
                          p3_electr_binary = mean(bg_sur.s3.l$p3_electr_binary, na.rm = T),
                          
                          temp_1 = mean(bg_sur.s3.l$temp_1, na.rm = T),
                          temp_2 = mean(bg_sur.s3.l$temp_2, na.rm = T),
                          temp_3 = mean(bg_sur.s3.l$temp_3, na.rm = T),
                          temp_4 = mean(bg_sur.s3.l$temp_4, na.rm = T),
                          temp_5 = mean(bg_sur.s3.l$temp_5, na.rm = T),
                          temp_6 = mean(bg_sur.s3.l$temp_6, na.rm = T),
                          temp_7 = mean(bg_sur.s3.l$temp_7, na.rm = T),
                          
                          education = "Baccalaureate",
                          date_month = median(bg_sur.s3.l$date_month),
                          date_week = median(bg_sur.s3.l$date_week),
                          mahala = "Kadhimyah",
                          block = "1001")

temperature_deviation_bg_logit <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$max > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + 
                                  temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                                  p3_electr_binary + 
                                  as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                                data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.l, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.l, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.l, new_data_bg[2,], type = "response") - 
                                    predict(reg_temp_week.l, new_data_bg[1,], type = "response")) /  mean(bg_sur.s3.l$vs_mnf)
    
    reg_temp_week_pars.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd +
                                       p3_electr_binary + 
                                       as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                                     data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.l, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.l, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + 
                                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                                   p3_electr_binary + 
                                   as.factor(education) + as.factor(block) + as.factor(date_month), 
                                 data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.l, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.l, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_month_pars.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd +
                                        p3_electr_binary + 
                                        as.factor(education) + as.factor(block) + as.factor(date_month), 
                                      data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$max)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.l, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.l, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$max)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

bg_95_logit <- temperature_deviation_bg_logit(bg_sur.s3.l, 0.95)
bg_90_logit <- temperature_deviation_bg_logit(bg_sur.s3.l, 0.90)

temp <- scales::percent(bg_95_logit[bg_95_logit$time %in% "week", ]$change, accuracy = 0.01)
temp <- replace(temp, seq(2,length(temp),2), "")
p_bg_logit <- ggplot() + 
  
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 1, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = bg_95_logit[bg_95_logit$time %in% "week", ]$gamma, labels=temp)) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0125, 0.055) + 
  annotate(geom = "text", x = 75:105, y = -0.0082,
           label = temp_ranges_bg,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  ggtitle("") +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"),
        plot.title=element_text(family='', hjust=0.0275, vjust=-6.75, face='bold', colour='black', size=14, margin=margin(t=40,b=-30))) +
  annotate(geom="text", x=75, y=0.05, label="Support for Violence",
           color="black", size = 6, fontface = 2, hjust = 0) + 
  annotate("point", x = 97, y = 0.049, colour = "blue4", size = 3) + 
  
  annotate(geom="text", x=98, y= 0.049, label="Primary Specification",
           color="black", size = 3, hjust = 0) + 
  annotate("point", x = 97, y = 0.0445, colour = "gray60", size = 3) + 
  annotate(geom="text", x=98, y= 0.0445, label="Robustness Checks",
           color="black", size = 3, hjust = 0) +
  geom_rect(aes(xmin = 95.5, xmax = 103, ymin = 0.0395, ymax = 0.0545),
            fill = "transparent", color = "black", size = 0.5)

# Plot:

baghdad_survey_combined <- ggarrange(p_bg_logit + xlab("") + ylab("Logistic Regression Model"),
                                     p_bg_ols + xlab("") + ylab("Linear Probability Model"),
                                     ncol = 1, nrow=2, labels  = "")

# PRODUCES FIGURE 23 IN PAPER:
annotate_figure(baghdad_survey_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Support Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

baghdad_survey_combined <- ggarrange(p_bg_logit + xlab("") + ylab("Logistic Regression Model"),
                                     ncol = 1, nrow=1, labels  = "")

# PRODUCES FIGURE 6 IN PAPER:
annotate_figure(baghdad_survey_combined,
                top = text_grob(TeX("$\\rho$: Estimated % Change in Support Following 1$\\sigma$ Increase in Temperature Deviation"), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*": High Temperature Threshold (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### END.
###

# Next, generate primary relationship using alternative approaches for testing for 
# non-linearity in the relationship. 

# Temperature bins:

# Iraq:
# Construct temperature bins (three-degree increments): 
iq_sig1$bin <- cut(iq_sig1$max, seq(min(round(iq_sig1$max)), max(round(iq_sig1$max)), by = 3))
# Then, regress the bins on violence:
IQ.w.u.ols.bin <- felm(log(unconstrained+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                         temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                         constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                         unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                         idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                         ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                         week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.ols.bin.ab <- felm(log(unconstrained+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                            week + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.ols.bin <- felm(log(unconstrained+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                         temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                         constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                         unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                         idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                         ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                         month + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.ols.bin.ab <- felm(log(unconstrained+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                            month + governorate | 0 | 0, data = iq_sig1)

stargazer(IQ.w.u.ols.bin, IQ.w.u.ols.bin.ab, IQ.m.u.ols.bin, IQ.m.u.ols.bin.ab,
          se = list(tidy(IQ.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error, 
                    tidy(IQ.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error,
                    tidy(IQ.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error, 
                    tidy(IQ.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error), type = "text")

df1.a <- tidy(IQ.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.90)
df1.b <- tidy(IQ.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)

df2.a <- tidy(IQ.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.90)
df2.b <- tidy(IQ.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)

df3.a <- tidy(IQ.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.90)
df3.b <- tidy(IQ.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)

df4.a <- tidy(IQ.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.90)
df4.b <- tidy(IQ.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)

bins_labels <- as.character(unique(sort(cut(iq_sig1$max, seq(min(round(iq_sig1$max)), max(round(iq_sig1$max)), by = 3)))))[2:26]

p.bins_iq <- ggplot() + 
  geom_point(data = df1.a[1:25, ], aes(y = estimate, x = 1:25), color = "blue4", size = 2) + 
  geom_errorbar(data = df1.a[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "blue4", linetype = "solid") +
  geom_errorbar(data = df1.b[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "blue4", linetype = "dashed") +
  
  geom_point(data = df2.a[1:25, ], aes(y = estimate, x = 1:25), color = "gray60") + 
  geom_errorbar(data = df2.a[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df2.b[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  geom_point(data = df3.a[1:25, ], aes(y = estimate, x = 1:25), color = "gray60") + 
  geom_errorbar(data = df3.a[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df3.b[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  geom_point(data = df4.a[1:25, ], aes(y = estimate, x = 1:25), color = "gray60") + 
  geom_errorbar(data = df4.a[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df4.b[1:25, ], aes(y = estimate, x = 1:25, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  scale_x_continuous(breaks = c(1:25), labels=c(bins_labels)) +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text.x = element_text(angle = 45, vjust=0.0)) +
  ylab(expression(beta)) + 
  xlab("")

# Afghanistan:
# Construct temperature bins (three-degree increments): 
af_sig1$bin <- cut(af_sig1$max, seq(min(round(af_sig1$max)), max(round(af_sig1$max)), by = 3))
# Then, regress the bins on violence:
AF.w.u.ols.bin <- felm(log(df+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                         temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                         df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                         idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                         ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                         week + district | 0 | 0, data = af_sig1)
AF.w.u.ols.bin.ab <- felm(log(df+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                            week + district | 0 | 0, data = af_sig1)
AF.m.u.ols.bin <- felm(log(df+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                         temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                         df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                         idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                         ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                         month + district | 0 | 0, data = af_sig1)
AF.m.u.ols.bin.ab <- felm(log(df+1) ~ as.factor(bin) + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                            month + district | 0 | 0, data = af_sig1)

stargazer(AF.w.u.ols.bin, AF.w.u.ols.bin.ab, AF.m.u.ols.bin, AF.m.u.ols.bin.ab,
          se = list(tidy(AF.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error, 
                    tidy(AF.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error,
                    tidy(AF.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error, 
                    tidy(AF.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)$std.error), type = "text")

df1.a <- tidy(AF.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.90)
df1.b <- tidy(AF.w.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)

df2.a <- tidy(AF.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.90)
df2.b <- tidy(AF.w.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)

df3.a <- tidy(AF.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.90)
df3.b <- tidy(AF.m.u.ols.bin, conf.int = T, se.type = "robust", conf.level = 0.95)

df4.a <- tidy(AF.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.90)
df4.b <- tidy(AF.m.u.ols.bin.ab, conf.int = T, se.type = "robust", conf.level = 0.95)

bins_labels <- as.character(unique(sort(cut(af_sig1$max, seq(min(round(af_sig1$max)), max(round(af_sig1$max)), by = 3)))))[4:35]

p.bins_af <- ggplot() + 
  geom_point(data = df1.a[1:32, ], aes(y = estimate, x = 1:32), color = "blue4", size = 2) + 
  geom_errorbar(data = df1.a[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "blue4", linetype = "solid") +
  geom_errorbar(data = df1.b[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "blue4", linetype = "dashed") +
  
  geom_point(data = df2.a[1:32, ], aes(y = estimate, x = 1:32), color = "gray60") + 
  geom_errorbar(data = df2.a[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df2.b[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  geom_point(data = df3.a[1:32, ], aes(y = estimate, x = 1:32), color = "gray60") + 
  geom_errorbar(data = df3.a[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df3.b[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  geom_point(data = df4.a[1:32, ], aes(y = estimate, x = 1:32), color = "gray60") + 
  geom_errorbar(data = df4.a[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  geom_errorbar(data = df4.b[1:32, ], aes(y = estimate, x = 1:32, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  
  scale_x_continuous(breaks = c(1:32), labels=c(bins_labels)) +
  theme_bw() +
  theme(text = element_text(size=12),
        axis.text.x = element_text(angle = 45, vjust=0.0)) +
  ylab(expression(beta)) + 
  xlab("")

# Baghdad Surveys:
bg_sur.s3.l$bin <- cut(bg_sur.s3.l$max, seq(min(round(bg_sur.s3.l$max)), max(round(bg_sur.s3.l$max)), by = 3))

BG.mm.u.lg.bin <- bayesglm(vs_mnf ~ as.factor(bin) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                             p3_electr_binary +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"), drop.unused.levels = F)
BG.mm.u.lg.bin.ab <- bayesglm(vs_mnf ~ as.factor(bin) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                                p3_electr_binary +
                                as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"), drop.unused.levels = F)
BG.mb.u.lg.bin <- bayesglm(vs_mnf ~ as.factor(bin) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                             p3_electr_binary +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"), drop.unused.levels = F)
BG.mb.u.lg.bin.ab <- bayesglm(vs_mnf ~ as.factor(bin) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                                p3_electr_binary +
                                as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"), drop.unused.levels = F)

# Then, construct a function for generating predicted probabilities and creating output that can be plotted.
plot_function <- function(reg_output, dataset, model){
  
  dataset <- dataset[!(is.na(dataset$max)), ]
  
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(reg_output), Sigma = vcov(reg_output))
  paste <- model.matrix(reg_output)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps.df1 <- matrix(NA, length(unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3))))))-1, 100)
  for(i in 1:length(unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3))))))-1){
    
    paste[, 2:length(unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3))))))] <- 0
    paste[, i+1] <- 1
    x.beta <- paste %*% t(beta.sim)
    
    yhat <- exp(x.beta) / (exp(x.beta) + 1)
    
    reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  plot_dataset <- as.data.frame(matrix(NA, length(unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3))))))-1,3))
  colnames(plot_dataset) <- c("means", "cil", "ciu")
  
  plot_dataset$means <- apply(reps.df1, 1 , mean)
  se <- apply(reps.df1, 1 , sd)
  plot_dataset$cil <- plot_dataset$means - qnorm(0.975)*se
  plot_dataset$ciu <- plot_dataset$means + qnorm(0.975)*se
  plot_dataset$temperature <- unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3)))))[2:length(unique(na.omit(sort(cut(dataset$max, seq(min(round(dataset$max)), max(round(dataset$max)), by = 3))))))]
  
  return(plot_dataset)
}

p.BG.mm.u.lg.bin <- plot_function(BG.mm.u.lg.bin, bg_sur.s3.l, "logit")
p.BG.mm.u.lg.bin.ab <- plot_function(BG.mm.u.lg.bin.ab, bg_sur.s3.l, "logit")
p.BG.mb.u.lg.bin <- plot_function(BG.mb.u.lg.bin, bg_sur.s3.l, "logit")
p.BG.mb.u.lg.bin.ab <- plot_function(BG.mb.u.lg.bin.ab, bg_sur.s3.l, "logit")

p.bins_bg <- ggplot() +
  geom_point(data = p.BG.mm.u.lg.bin, aes(y=means, x = temperature)) + 
  geom_errorbar(data = p.BG.mm.u.lg.bin, aes(y = means, x = temperature, ymin = cil, ymax = ciu), alpha = .35, size = 0.4, color = "blue4", linetype = "solid") +
  
  geom_point(data = p.BG.mm.u.lg.bin.ab, aes(y=means, x = temperature), color = "gray60") + 
  geom_errorbar(data = p.BG.mm.u.lg.bin.ab, aes(y = means, x = temperature, ymin = cil, ymax = ciu), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  
  geom_point(data = p.BG.mb.u.lg.bin, aes(y=means, x = temperature), color = "gray60") + 
  geom_errorbar(data = p.BG.mb.u.lg.bin, aes(y = means, x = temperature, ymin = cil, ymax = ciu), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") +
  
  geom_point(data = p.BG.mb.u.lg.bin.ab, aes(y=means, x = temperature), color = "gray60") + 
  geom_errorbar(data = p.BG.mb.u.lg.bin.ab, aes(y = means, x = temperature, ymin = cil, ymax = ciu), alpha = .35, size = 0.4, color = "gray60", linetype = "solid") + 
  
  theme_bw() + 
  
  theme(text = element_text(size=12),
        axis.text.x = element_text(angle = 45, vjust=0.0)) +
  ylab("Predicted Probability") + 
  xlab("")

# Combine the figures:
bins_combined <- ggarrange(p.bins_iq, p.bins_af, p.bins_bg, ncol = 3)

bins_combined <- annotate_figure(bins_combined,
                                 left = text_grob("Binned Regression", color = "Black", rot = 90, size = 16),
                                 bottom = text_grob(expression("Bins (°F)"), color = "Black", rot = 0, size = 16))

###
###
###

# General Additive Models:

require(devtools)
install_version("mgcv", version = "1.8.1", repos = "http://cran.us.r-project.org")

# Iraq:

IQ.w.u.gam <- gam(unconstrained ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + hoursofpower + idf + ied + constrained +
                    temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                    constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                    unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                    idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                    ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                    as.factor(week) + as.factor(governorate), 
                  data = iq_sig1, family = poisson)

IQ.w.u.gam.ab <- gam(unconstrained ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + hoursofpower + idf + ied + constrained +
                       as.factor(week) + as.factor(governorate), 
                     data = iq_sig1, family = poisson)

IQ.m.u.gam <- gam(unconstrained ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + hoursofpower + idf + ied + constrained +
                    temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                    constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                    unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                    idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                    ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                    as.factor(month) + as.factor(governorate), 
                  data = iq_sig1, family = poisson)

IQ.m.u.gam.ab <- gam(unconstrained ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + hoursofpower + idf + ied + constrained +
                       as.factor(month) + as.factor(governorate), 
                     data = iq_sig1, family = poisson)

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(IQ.w.u.gam)
1
1
1
1
1

IQ.results.gam <- as.data.frame(plotData[[1]]$fit)
colnames(IQ.results.gam) <- "r1"
IQ.results.gam$max <- plotData[[1]]$x
IQ.results.gam$upper1 <- plotData[[1]]$fit + plotData[[1]]$se
IQ.results.gam$lower1 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(IQ.w.u.gam.ab)
1
1
1
1
1

IQ.results.gam$r2 <- plotData[[1]]$fit
IQ.results.gam$upper2 <- plotData[[1]]$fit + plotData[[1]]$se
IQ.results.gam$lower2 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(IQ.m.u.gam)
1
1
1
1
1

IQ.results.gam$r3 <- plotData[[1]]$fit
IQ.results.gam$upper3 <- plotData[[1]]$fit + plotData[[1]]$se
IQ.results.gam$lower3 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(IQ.m.u.gam.ab)
1
1
1
1
1

IQ.results.gam$r4 <- plotData[[1]]$fit
IQ.results.gam$upper4 <- plotData[[1]]$fit + plotData[[1]]$se
IQ.results.gam$lower4 <- plotData[[1]]$fit - plotData[[1]]$se

p.gam.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam, aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam, aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam, aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam, aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam, aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam, aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam, aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam, aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  scale_x_continuous(breaks=seq(min(round(IQ.results.gam$max)), max(round(IQ.results.gam$max)), by = 5)) + 
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Small Arms Fire, Iraq")

# Afghanistan:

AF.w.u.gam <- gam(df ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + idf + ied +
                    temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                    df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                    idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                    ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                    as.factor(week) + as.factor(district), 
                  data = af_sig1, 
                  family = poisson)
AF.w.u.gam.ab <- gam(df ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + idf + ied +
                       as.factor(week) + as.factor(district), 
                     data = af_sig1, 
                     family = poisson)
AF.m.u.gam <- gam(df ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + idf + ied +
                    temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                    df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                    idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                    ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                    as.factor(month) + as.factor(district), 
                  data = af_sig1, 
                  family = poisson)
AF.m.u.gam.ab <- gam(df ~ s(max) + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight + idf + ied +
                       as.factor(month) + as.factor(district), 
                     data = af_sig1, 
                     family = poisson)

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(AF.w.u.gam)
1
1
1
1
1

AF.results.gam <- as.data.frame(plotData[[1]]$fit)
colnames(AF.results.gam) <- "r1"
AF.results.gam$max <- plotData[[1]]$x
AF.results.gam$upper1 <- plotData[[1]]$fit + plotData[[1]]$se
AF.results.gam$lower1 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(AF.w.u.gam.ab)
1
1
1
1
1

AF.results.gam$r2 <- plotData[[1]]$fit
AF.results.gam$upper2 <- plotData[[1]]$fit + plotData[[1]]$se
AF.results.gam$lower2 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(AF.m.u.gam)
1
1
1
1
1

AF.results.gam$r3 <- plotData[[1]]$fit
AF.results.gam$upper3 <- plotData[[1]]$fit + plotData[[1]]$se
AF.results.gam$lower3 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(AF.m.u.gam.ab)
1
1
1
1
1

AF.results.gam$r4 <- plotData[[1]]$fit
AF.results.gam$upper4 <- plotData[[1]]$fit + plotData[[1]]$se
AF.results.gam$lower4 <- plotData[[1]]$fit - plotData[[1]]$se

# For presentation purposes, restrict lower confidence interval to -5:
AF.results.gam$lower2 <- ifelse(AF.results.gam$lower2 < -5, -5, AF.results.gam$lower2)

p.gam.AF <- ggplot() + 
  geom_line(data = AF.results.gam, aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam, aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam, aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam, aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam, aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam, aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam, aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam, aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  scale_x_continuous(breaks=seq(min(round(AF.results.gam$max)), max(round(AF.results.gam$max)), by = 5)) + 
  
  xlab("") + 
  ylab("") +
  ylim(-7,2) +
  ggtitle("Direct Fire, Afghanistan")

# Baghdad surveys:

BG.mm.u.gam <- gam(vs_mnf ~ s(max) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight +
                     p3_electr_binary +
                     temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                     as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                   data = bg_sur.s3.l, 
                   family=binomial, method = "REML")
BG.mm.u.gam.ab <- gam(vs_mnf ~ s(max) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight +
                        p3_electr_binary +
                        as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                      data = bg_sur.s3.l, 
                      family=binomial, method = "REML")
BG.mb.u.gam <- gam(vs_mnf ~ s(max) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight +
                     p3_electr_binary +
                     temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                     as.factor(education) + as.factor(block) + as.factor(date_month), 
                   data = bg_sur.s3.l, 
                   family=binomial, method = "REML")
BG.mb.u.gam.ab <- gam(vs_mnf ~ s(max) + age + work_wkhr + income.w + hhld_size + m_sig1_7 + s(prcp) + s(wdsp) + s(dewp) + s(visib) + s(mxspd) + daylight +
                        p3_electr_binary +
                        as.factor(education) + as.factor(block) + as.factor(date_month), 
                      data = bg_sur.s3.l, 
                      family=binomial, method = "REML")

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(BG.mm.u.gam)
1
1
1
1
1

BG.results.gam <- as.data.frame(plotData[[1]]$fit)
colnames(BG.results.gam) <- "r1"
BG.results.gam$max <- plotData[[1]]$x
BG.results.gam$upper1 <- plotData[[1]]$fit + plotData[[1]]$se
BG.results.gam$lower1 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(BG.mm.u.gam.ab)
1
1
1
1
1

BG.results.gam$r2 <- plotData[[1]]$fit
BG.results.gam$upper2 <- plotData[[1]]$fit + plotData[[1]]$se
BG.results.gam$lower2 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(BG.mb.u.gam)
1
1
1
1
1

BG.results.gam$r3 <- plotData[[1]]$fit
BG.results.gam$upper3 <- plotData[[1]]$fit + plotData[[1]]$se
BG.results.gam$lower3 <- plotData[[1]]$fit - plotData[[1]]$se

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
      quote({plotData <<- pd}))
mgcv::plot.gam(BG.mb.u.gam.ab)
1
1
1
1
1

BG.results.gam$r4 <- plotData[[1]]$fit
BG.results.gam$upper4 <- plotData[[1]]$fit + plotData[[1]]$se
BG.results.gam$lower4 <- plotData[[1]]$fit - plotData[[1]]$se

p.gam.BG <- ggplot() + 
  geom_line(data = BG.results.gam, aes(x = max, y = r1)) + 
  geom_ribbon(data =  BG.results.gam, aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = BG.results.gam, aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  BG.results.gam, aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = BG.results.gam, aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  BG.results.gam, aes(x = max, y = r3, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = BG.results.gam, aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  BG.results.gam, aes(x = max, y = r4, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  scale_x_continuous(breaks=seq(min(round(BG.results.gam$max)), max(round(BG.results.gam$max)), by = 5)) + 
  
  xlab("") + 
  ylab("") +
  ggtitle("Support for Attacks")

# Combine the plots:
gam_combined <- ggarrange(p.gam.IQ, p.gam.AF, p.gam.BG, ncol = 3)

gam_combined <- annotate_figure(gam_combined,
                                left = text_grob("Generalized Additive Models", color = "Black", rot = 90, size = 16),
                                bottom = text_grob(expression("Temperature (°F)"), color = "Black", rot = 0, size = 16))

# PRODUCES FIGURE 19 IN PAPER:
ggarrange(gam_combined, bins_combined, nrow = 2)

###
###
###

# Next, show alternative weather associations:
# Iraq:
IQ.results.gam_alt <- list()
IQ.results.gam_alt[[1]] <- NA

for(i in 2:6){
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(IQ.w.u.gam)
  
  
  IQ.results.gam_alt[[i]] <-  as.data.frame(plotData[[i]]$fit)
  colnames(IQ.results.gam_alt[[i]]) <- "r1"
  IQ.results.gam_alt[[i]]$max <- plotData[[i]]$x
  IQ.results.gam_alt[[i]]$upper1 <- plotData[[i]]$fit + plotData[[i]]$se
  IQ.results.gam_alt[[i]]$lower1 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(IQ.w.u.gam.ab)
  
  
  IQ.results.gam_alt[[i]]$r2 <- plotData[[i]]$fit
  IQ.results.gam_alt[[i]]$upper2 <- plotData[[i]]$fit + plotData[[i]]$se
  IQ.results.gam_alt[[i]]$lower2 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(IQ.m.u.gam)
  
  
  IQ.results.gam_alt[[i]]$r3 <- plotData[[i]]$fit
  IQ.results.gam_alt[[i]]$upper3 <- plotData[[i]]$fit + plotData[[i]]$se
  IQ.results.gam_alt[[i]]$lower3 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(IQ.m.u.gam.ab)
  
  
  IQ.results.gam_alt[[i]]$r4 <- plotData[[i]]$fit
  IQ.results.gam_alt[[i]]$upper4 <- plotData[[i]]$fit + plotData[[i]]$se
  IQ.results.gam_alt[[i]]$lower4 <- plotData[[i]]$fit - plotData[[i]]$se
}
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

p.prcp.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam_alt[[2]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam_alt[[2]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam_alt[[2]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[2]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[2]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[2]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[2]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[2]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Precipitation")

p.wdsp.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam_alt[[3]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam_alt[[3]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam_alt[[3]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[3]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[3]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[3]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[3]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[3]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Wind Speed")

p.dewp.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam_alt[[4]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam_alt[[4]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam_alt[[4]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[4]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[4]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[4]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[4]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[4]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Dew Point")

p.visib.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam_alt[[5]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam_alt[[5]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam_alt[[5]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[5]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[5]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[5]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[5]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[5]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Visibility")

p.mxspd.IQ <- ggplot() + 
  geom_line(data = IQ.results.gam_alt[[6]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  IQ.results.gam_alt[[6]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = IQ.results.gam_alt[[6]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[6]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[6]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[6]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = IQ.results.gam_alt[[6]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  IQ.results.gam_alt[[6]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Maximum Wind Speed")

##
##

# Afghanistan:

AF.results.gam_alt <- list()
AF.results.gam_alt[[1]] <- NA

for(i in 2:6){
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(AF.w.u.gam)
  
  
  AF.results.gam_alt[[i]] <-  as.data.frame(plotData[[i]]$fit)
  colnames(AF.results.gam_alt[[i]]) <- "r1"
  AF.results.gam_alt[[i]]$max <- plotData[[i]]$x
  AF.results.gam_alt[[i]]$upper1 <- plotData[[i]]$fit + plotData[[i]]$se
  AF.results.gam_alt[[i]]$lower1 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(AF.w.u.gam.ab)
  
  
  AF.results.gam_alt[[i]]$r2 <- plotData[[i]]$fit
  AF.results.gam_alt[[i]]$upper2 <- plotData[[i]]$fit + plotData[[i]]$se
  AF.results.gam_alt[[i]]$lower2 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(AF.m.u.gam)
  
  
  AF.results.gam_alt[[i]]$r3 <- plotData[[i]]$fit
  AF.results.gam_alt[[i]]$upper3 <- plotData[[i]]$fit + plotData[[i]]$se
  AF.results.gam_alt[[i]]$lower3 <- plotData[[i]]$fit - plotData[[i]]$se
  
  plotData <- list()
  trace(mgcv:::plot.gam, at = list(c(27, 1)), 
        quote({plotData <<- pd}))
  mgcv::plot.gam(AF.m.u.gam.ab)
  
  
  AF.results.gam_alt[[i]]$r4 <- plotData[[i]]$fit
  AF.results.gam_alt[[i]]$upper4 <- plotData[[i]]$fit + plotData[[i]]$se
  AF.results.gam_alt[[i]]$lower4 <- plotData[[i]]$fit - plotData[[i]]$se
}
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

p.prcp.AF <- ggplot() + 
  geom_line(data = AF.results.gam_alt[[2]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam_alt[[2]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam_alt[[2]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[2]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[2]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[2]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[2]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[2]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Precipitation")

p.wdsp.AF <- ggplot() + 
  geom_line(data = AF.results.gam_alt[[3]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam_alt[[3]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam_alt[[3]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[3]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[3]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[3]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[3]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[3]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Wind Speed")

p.dewp.AF <- ggplot() + 
  geom_line(data = AF.results.gam_alt[[4]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam_alt[[4]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam_alt[[4]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[4]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[4]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[4]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[4]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[4]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Dew Point")

p.visib.AF <- ggplot() + 
  geom_line(data = AF.results.gam_alt[[5]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam_alt[[5]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam_alt[[5]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[5]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[5]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[5]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[5]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[5]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Visibility")

p.mxspd.AF <- ggplot() + 
  geom_line(data = AF.results.gam_alt[[6]], aes(x = max, y = r1)) + 
  geom_ribbon(data =  AF.results.gam_alt[[6]], aes(x = max, y = r1, ymin = lower1, ymax = upper1), alpha = .15, fill = "blue4") +
  
  geom_line(data = AF.results.gam_alt[[6]], aes(x = max, y = r2), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[6]], aes(x = max, y = r2, ymin = lower2, ymax = upper2), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[6]], aes(x = max, y = r3), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[6]], aes(x = max, y = r3, ymin = lower3, ymax = upper3), alpha = .15, fill = "blue4") + 
  
  geom_line(data = AF.results.gam_alt[[6]], aes(x = max, y = r4), color = "gray60") + 
  geom_ribbon(data =  AF.results.gam_alt[[6]], aes(x = max, y = r4, ymin = lower4, ymax = upper4), alpha = .15, fill = "blue4") + 
  theme_bw() +
  
  theme(axis.ticks = element_blank(), 
        text = element_text(size=12),
        axis.text.x = element_text(angle = 45)) +
  
  xlab("") + 
  ylab(expression(beta)) +
  ggtitle("Maximum Wind Speed")

# Export figures:

alt_IQ <- ggarrange(p.prcp.IQ, p.wdsp.IQ, p.mxspd.IQ, p.dewp.IQ, p.visib.IQ, nrow = 5)
alt_IQ <- annotate_figure(alt_IQ,
                          top = text_grob("Iraq", color = "Black", rot = 0, size = 16))

alt_AF <- ggarrange(p.prcp.AF, p.wdsp.AF, p.mxspd.AF, p.dewp.AF, p.visib.AF, nrow = 5)
alt_AF <- annotate_figure(alt_AF,
                          top = text_grob("Afghanistan", color = "Black", rot = 0, size = 16))

# PRODUCES FIGURE 18 IN PAPER:
ggarrange(alt_IQ, alt_AF, ncol = 2)

###
### END.
###
